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ABSTRACT 


This  report  summarizes  the  results  of  research  during 
the  past  two  years,  in  understanding  the  resolution  of  the 
autoregressive  (AR)  spectral  estimators  and  developing  and 
evaluating  computationally  efficient  autoregressive-moving 
average  (ARMA)  spectral  estimators.  The  loss  in  the  resolution 
of  the  AR  spectral  estimator  in  the  presence  of  noise  is  related 
to  the  appearance  of  zeros  in  the  z-plane.  A  parallel  resonator 
model  is  developed  to  relate  the  loss  in  resolution  (bandwidth 
expansion)  to  the  signal- to-noise  ratio  and  parameters  of  the 
noiseless  signal  model. 

A  new  technique  for  the  identification  of  the  order  of 
an  AR  model  was  derived  that  shows  substantial  stability  compared 
to  the  popular  Akaike  Information  Criterion  method.  Order 
determination  was  emphasized,  since  increase  in  the  order  of 
the  AR  spectral  estimator,  to  account  for  the  presence  of  noise, 
is  naturally  accompanied  by  larger  variance  of  the  estimates  and 
appearance  of  spurious  peaks. 

Several  sub-optimum  (non-maximum-likelihood)  ARMA  spectral 
estimators  were  also  developed.  These  methods  are  computationally 
efficient,  but  statistically  not  very  stable  for  small  data  records. 
An  evaluation  of  the  statistical  properties  of  the  different  sub¬ 
optimum  ARMA  techniques  led  to  the  evaluation  of  asymptotic  bounds 
on  the  variances  of  the  estimates  of  the  parameters  or  the  poles 
and  zeros  of  the  model  through  the  evaluation  of  Fisher's 
information  matrix.  Finally,  a  modification  of  Burg's  MEM  spectral 
estimator  was  developed  that  improves  the  accuracy  of  the  spectral 
estimates  of  complex  sinusoids  and  makes  the  method  considerably 


more  robust. 
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I .  INTRODUCTION 

The  problem  of  spectral  estimation  via  the  autoregressive 
(AR)  and  autoregressive-moving-average  (ARMA)  modeling  of  the 
observed  noisy  signal  is  treated  in  this  report.  High  resolution 
spectral  analysis  based  on  the  AR  model  has  received  much 
attention  in  recent  years  [  1  ]  -  [  3  ]  .  The  basic  impetus  for 

taking  these  approaches  is  due  to  the  fact  that  the  parametric 
modeling  schemes  are  data  adaptive  and  are  free  from  the  effects 
of  window  functions  that  a’'.:  inherent  in  the  traditional 
Blackman  and  Tukey  [  4  ]  type  spectral  estimators.  Furthermore, 
the  AR  model  is  easily  estimated  making  it  useful  in  applications 
such  as  radar  signal  processing  that  require  near- real- time 
processing . 

The  properties  of  the  AR  spectral  estimator  have  been 
studied,  theoretically  in  the  asymptotic  case  [  3  ]  ,  [  5  ]  ,  [  6  ] 
and  empirically  [  1  ]  ,  [  2  } .  It  has  been  shown  that  this 
estimator  in  many  cases  offers  considerably  higher  resolution 
based  on  the  same  amount  of  data,  than  the  Blackman  and  Tukey 
type  estimators.  Furthermore,  the  bove  asymptotic  and  empirical 
investigations  have  shown  the  variance  of  the  AR  estimates  to 
be  comparable  to  the  unsmoothed  Blackman  and  Tukey  type  estimates, 
for  the  same  number  of  autocorrelation  lags.  It  should  be 
pointed  out,  however,  that  the  AR  estimates  usually  reouire  much 
fewer  lags  for  the  same  resolution. 

Most  practical  applications  of  spectral  analysis  involve 
noisy  signals  and/or  multiple  signal  and  noise  mixtures. 

Therefore,  it  is  natural  to  consider  the  performance  of  the  AR 
spectral  estimator  in  the  presence  of  noise.  Previous  studies 
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[13/  (  2  ]  had  shown  resolution  degradation  in  the  presence 

of  noise.  The  solution  presented  for  improving  resolution  was 
given  as  increasing  the  order  of  the  estimator,  in  a  rather 
arbitrary  fashion. 

This  report  deals  with  the  question  of  AR  signals  in  noise. 

The  main  thrust  of  the  work  reported  here  was  to  improve  spectral 
resolution  by  using  ARMA  models  for  the  measured  signals. 
Specifically  some  suboptimum  but  computationally  efficient  ARMA 
spectral  estimation  algorithms  were  developed  and  their  properties 
studied.  In  the  process  of  the  investigations,  an  order 
determination  scheme  as  well  as  a  new  Burg  type  spectral 
estimator,  however,  were  also  developed  and  will  be  described  here, 

1.1  Report  Outline 

This  report  is  organized  as  follows.  First  the  AR  process 
is  defined.  The  sum  of  uncorrelated  AR  processes  and  white 


noise  is  then  considered  and  shown  to  be  represented  by  an 
equivalent  ARMA  process.  A  multiple  resonator  model  for  signals 
of  interest  is  presented  and  shown  to  be  equivalent  to  the  ARMA 
model.  Resolution  degradation  as  a  function  of  noise  and 
resonator  pole  locations  is  considered  and  from  it  some 
representative  curves  of  bandwidth  expansion  are  presented. 

Chapter  III.  of  the  report  deals  with  the  basic  question 
of  AR  model  order  determination.  A  new  criterion  related  to  that 
of  Akaike's  (AIC)  is  derived  and  comparative  examples  are  given. 
The  next  two  chapters  treat  the  basic  emphasis  of  this  work, 
namely  that  of  ARMA  spectral  estimation.  First,  several  ARMA 
schemes  are  presented  that  are  computationally  efficient  but  sub¬ 
optimum.  Sub-optimality  is  treated  with  respect  to  the  maximum- 
likelihood  (ML)  estimates  of  the  parameters.  The  subsequent 


chapter  then  discusses  some  statistical  properties  of  classes 
of  ARMA  spectral-estimator,  with  respect  to  the  statistic  used 
in  the  estimation.  Here,  again,  the  reference  scheme  is  the 
ML  estimator.  The  final  chapter  of  the  report  deals  with  a 
new  scheme,  related  to  Burg's  MEM  method,  which  is  especially 
useful  for  the  common  case  of  sinusoidal  signals  in  noise.  This 
method  is  shown  to  be  more  robust  than  the  MEM  technique,  with 
better  accuracy  and  equivalent  or  superior  resolution. 


I  iiu  t  i  iVs 
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II.  AUTOREGRESSIVE  SPECTRAL 
ESTIMATION  OF  NOISY  SIGNALS 

The  most  popular  recent  data-adaptive  spectral  estimation 
method  is  one  based  on  an  all-pole  model  for  the  siqnal.  The 
algorithmic  approach  for  the  estimation  of  the  parameters  of 
such  a  system  include,  the  methods  of  fitting  the  autoregressive 
(AR)  coefficients  to  the  data  as  well  as  the  popular  Burg  maximum 
entropy  method  MEM.  Because  of  the  popularity  of  this  model 
and  computational  simplicity  of  its  estimates,  the  noisy  signal 
spectral  estimation  will  be  confined  to  autoregressive  signal 
models.  That  is,  it  will  be  assumed  that  in  the  absence  of 
noise  and  interference  an  AR  model  satisfactorily  describes 
the  signal.  Resolution  deqradation  of  the  AR  spectral  estimates 
in  the  presence  of  noise  and  interference  is  then  investigated 
and  related  to  the  changes  in  the  model  structure. 

II. 1  The  AR  Spectral  Estimator 

+■  Vi 

A  zero-mean  time  series  {x*.}  is  said  to  satisfy  an  l  order 
autoregressive  model  if: 

L 

Xt  =  ai  xt-i  +  ut 

i  =  l 


where  {a^}  denote  the  AR  coefficients  and  {ut}  is  a  zero-mean 
uncorrelated  (white)  sequence.  Another  interpretation  of  the 
model  in  (II.  1)  is  that  {a^}  represent  an  Lt^'  order  one-step  ahead 
linear  predictor  of  {x^}.  If  are  then  estimated,  based  on 

a  minimum  mean  square  error  criterion  {ufc}  on  the  average  becomes 
an  orthogonal  sequence.  It  can  be  shown  [2  ]  that  the  model  in 
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^  i2irkf  AT  i  2 

E  a  ■  e  | 

k=l 


(II. 2) 


where  S-^  is  the  spectral  level  of  the  {utl  sequence. 

The  spectral  density  shown  in  (II. 2)  is  the  AR  spectrum 
and  it  is  this  model  that  is  fitted  to  an  observed  time  series, 
by  simply  estimating  {a^}  from  the  time  series.  Several  estimation 
procedures  for  {a^}  have  been  discussed  in  the  literature  such 
as  the  maximum  likelihood,  the  least-squares  [7]  and  Burg's 
method  based  on  forward  and  backward  prediction  error  filtering 
[8].  If  the  number  of  data  samples  is  not  very  small, 
the  simplest  and  computational ly  most  efficient  estimates  of 
{a^}  are  the  solution  of  the  Yule-Walker  equations.  These 
equations  arise,  simply  by  multiplying  equation  (II. 1)  by 
xt_^,  i=l,...L  and  taking  the  expectation  of  both  sides,  to 
obtain  a  relation  between  the  autocorrelation  function  of  the 
process  {x._ }  and  the  coefficients  {a^}.  The  Yule-Walker  equations 
are  then  given  by: 

R  A  =  (II. 3) 

o  o 

where 


ro  rl  -  rL-l 

rl 

al 

ri 

R  = 

,  A  = 

and  = 

o 

rL- 1  ro 

• 

ql 

o 

rL 

t  h 

and  where  r^  is  the  autocorrelation  function  of  {xfc}  at  the  i  lag. 
Furthermore,  the  power  in  {ut}  can  be  found  by  multiplying  (II. 1) 


by  xt  and  taking  the  expectation  as: 


(II. 4) 


S,  =  [r  -  E  a-r.]AT  ,  AT  the  sampling  interval 
1  °  i=l  1  1 

In  practice  {cu}  and  are  estimated  from  (II.  3)  and  (II. 4)  based 
on  estimates  of  the  autocorrelation  function  {r^}. 

II. 2  Spectral  Estimation  in  the  Presence  of  Noise  and  Interference 
We  now  assume  that  the  signal  {xfc}  satisfies  an  L  h  order 
AR  model  and  therefore  its  spectrum  can  be  estimated  as  in  the 
previous  section.  The  problem  of  interest  is  the  estimation  of 
the  spectrum  of  the  observed  signal 


Yt  =  +  wt  +  nfc 


(II. 5) 


where  }  is  an  AR (M)  process  and  considered  to  be  the  inter¬ 
ference  and  (nfc}  is  a  white  noise  sequence,  with  {nfc},  {xfc}  and 
{ uj t }  mutually  uncorrelated. 

The  resolution  of  the  spectral  estimators  are  now  discussed  in 
the  asymptotic  case,  that  is,  when  the  autocorrelation  function 
of  {y  }  is  accurately  known. 

Let  { }  be  described  by  (II.  1)  and  fwt}  be  qiven  by  the 
following  autoregressive  model 
M 

u).  =  n  b .  ui.  .  +  v. 

t  .  i  t-i  t 

1=1 

The  z-spectrum  of  yfc  is  then  given  by: 


Sy(z)  = 


-1, 


D  (z) D  (z  )  D  (z)  D  (z  ) 

X  X  U)  a) 


—  +  Sn 


(II. 6) 


where  and  are  given  by  (II.  4)  using  the  appropriate  auto¬ 
correlation  values  for  {xfc}  and  {wt),  Sn  is  the  spectrum  of  nfc  and 
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L’l  . 

D  ( z )  =  1  -  E  a.z1  ,  D  (z)  =  1  -  Z  b.z1 
x  i=l  1  ^  i=l  1 

Putting  (II. 6)  under  a  common  denominator,  it  becomes  obvious 
that  Sy(z)  is  the  spectrum  of  an  autoregressive  moving  average 
process  of  orders  L+M  and  L+M,  i.e.,  AR ( L+M) /MA ( L+M) .  This  is 
a  process  with  L+M  zeros  and  L+M  poles,  where,  from  equation 
(II. 6),  the  numerator  polynomial  coefficients  are  related  to 
{a^},  {b^}  and  Sn  in  an  obvious  manner. 

It  can  be  seen  that  the  estimation  of  Sy(z)  using  a  purely 
AR  technique  (all-pole)  is  equivalent  to  approximating  the  (L+M) 
order  moving  average  component  by  an  AR  one.  This,  theoretically 
would  require  an  infinite  order  model.  Finite  order  models  of 
order  L,  however,  will  estimate  Sy(f)  with  good  resolution,  with 
L  depending  on  the  various  parameters  of  signal  noise  and  inter¬ 
ference,  notably  their  relative  power. 

It  is  obvious  now  that  whereas  the  resolution  of  Blackman 
and  Tukey  type  spectral  estimations  are  only  determined  by  the 
window  bandwidths  in  a  predictable  fashion,  those  of  the  AR  and 
generally  ARMA  estimators  are  very  much  data  dependent,  requirina 
larger  all-pole  orders  or  ARMA  modeling.  Figure  (II. 1)  shows  the 
spectrum  of  a  noiseless  AR  signal.  Figure  (II.  2)  shows  the 
calculated  AR  (from  exact  values  of  the  autocorrelation  function) 
spectrum  of  the  same  signal  in  the  presence  of  white  noise,  using 
order  L  =  20.  The  degradation  of  spectral  resolution  is  obvious. 
A  different  approach  to  the  demonstration  of  the  dependence  of 
spectral  resolution  on  the  signal- to-noise  ratio  is  to  consider  a 
parallel  resonator  model  for  the  measured  signal.  This  is 
discussed  in  the  next  section. 
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II. 3  Parallel  Resonator  Model 

High  resolution  spectral  estimation  is  normally  used  in 
situations  where  the  signal  has  a  "peaky"  spectrum.  Therefore, 
one  may  postulate  the  signal  model  as  the  sum  of  the  outputs  of 
M  second-order  resonators  driven  by  white  noise.  This  model, 
as  will  be  shown  in  the  sequel,  lends  itself  to  an  investigation 
of  the  resolution  degradation  or  bandwidth  expansion  as  a  function 
of  measurement  noise. 

Mathematically,  the  signal  s(n)  is  modeled  as  [  9  ], 

M 

s(n)  =  1  s  (n) 

m=l 

where 

sm(n)  =  +  am2sm(n'2)  +  am0u(n)  (II*7) 

where  u(n)  is  the  driving  noise  of  the  signal  generation  process. 

2 

It  is  assumed  that  u(n)  is  white  with  variance  oy.  The  signal 
model  is  shown  in  figure  (II. 3). 

The  transfer  function  of  the  signal  generation  process 
is  defined  as 


H  (z) 


S  ( z ) 

uTST 


(II. 8) 


where  S(z)  and  U(z)  denote  the  z-transforms  of  the  signal  s(n) 
and  the  noise  u(n)  respectively.  The  transfer  function  Hm(z)  of 
the  m'th  resonator  is  obtained  by  taking  the  z- transform  of  sm(n) 

Sm(z)  =  amlz"lsm(z)  +  am2z"2sm(z)  +  am0U(z) 


Hence , 


H  (z) 
m 


Sm(z) 

U(z) 


am0 


-  a 


ml 


-1 

z 


m2 


-2  * 

z 


1 


(II. 9) 


FIGURE  11.3  THE  SIGNAL  MODEL 
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The  z-transform  of  the  signal  s(n)  is  then 


S(z)  =  Z  Sm(z> 
m=l 


=  Z  U(z)H  (z) 

m=l 


We  can,  therefore,  express  the  transfer  function  for  the 


generation  of  s(n)  as: 


H(z)  = 


MM  12 

ak0  (1  ‘  amlZ"  _  am2z”  ) 

K=i  m=l 


- 1  -  2 
II  ( 1  ~  a  .  z  -  a  _  z  ) 
.  ml  m2 

m=l 


(II . 10) 


It  can  easily  be  seen  that  the  order  of  the  numerator  polynomial 
N(z)  of  the  transfer  function  H(z)  is  2M-2  while  the  order  of 
the  denominator  polynomial  is  2M.  H(z)  is  also  recognized  as 
the  transfer  function  of  an  autoregressive-moving  average  (ARMA) 
model.  Therefore,  using  the  parallel  resonator  model  is  the  same 
as  using  an  ARMA  model  where  the  difference  in  order  between  the 
autoregressive  (AR)  and  moving  average  (MA)  terms  is  two. 

Thus,  in  the  absence  of  noise,  s(n)  can  accurately  be  modeled 
as  an  ARMA  process  (or  the  output  of  a  pole-zero  filter) .  If  s(n) 
is  corrupted  by  the  uncorrelated  white  noise  sequence  w(n),  the 
z-power  spectrum  of  the  observation,  y(n),  is  given  by 

-  -  |N(z)  2  2  2 

y  ] D (z)  u  w 

cr2  |  N(z)  |  2  +  a 2 | D ( z ) | 2 


therefore,  in  the  presence  of  noise,  the  model  of  y(n)  is 
equivalent  to  an  ARMA  (2M,  2M)  process.  This  of  course  is 
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consistent  with  the  model  derived  in  the  preceding  section 
directly  from  the  addition  of  AR  processes  and  white  noise. 


II. 4  Analysis  of  a  Single  Resonator 

Some  insight  into  the  behavior  of  the  parallel  resonator 
model  in  the  presence  of  noise  can  be  gained  by  examining  the 
behavior  of  each  resonator  individually. 

To  simplify  the  analysis,  we  assume  that  the  noise  w(n)  is 
added  equally  to  the  output  of  each  resonator.  In  the  following 
first  the  set  of  conditions  which  the  parameters  of  each  resonator 
must  satisfy  for  stability  and  resonance  are  discussed.  The 
relations  between  the  parameters  and  the  location  of  the  poles 
of  a  resonator  are  then  derived.  The  effect  of  additive  noise 
on  the  poles  is  then  examined.  Finally,  the  perturbations  of 
the  pole  positions  of  the  noisy  resonator  as  a  function  of  the 
signal- to-noise  ratio  is  analyzed. 

The  transfer  function  of  the  m'th  resonator  is  given  by 

Sm(z) 

Hm<z)  =  UTFT 


mO 


-1  -2 
a  ,  z  -  a  z 
ml  m2 


(II. ID 


mO 


A(z) 

m 


The  poles  of  H  (z)  must  be  complex  conjugates  for  S  (n)  to  be  a 
m  m 

narrow  uand  process.  The  locations  of  the  poles  are  the  solution 
of  the  quadratic  equation 


2  2 

2  A  ( z )  =  z  -  a.z-ao  =  0 
m  ml  m2 


(11.12) 


or 
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z  = 


ml 


aml  +  4am2 


the  condition  for  having  complex  conjugate  poles  is  that 
the  discriminant  must  always  be  negative.  That  is. 


or 


+  4a 


m2 


<  0 


(11.13) 


If  equation  (11.13)  is  satisfied,  then  the  poles  would  be  in 
polar  coordinates, 

jo 

Z  —  Cj  0 

m  m 

where  om  is  the  distance  from  the  origin  of  the  unit  circle  to 
location  of  the  pole  and  tom  is  the  frequency  of  resonance  in 
radians/second. 

.  * 

The  system  is  stable  if  z  and  z  are  located  within  the 

m  ir. 

unit  circle.  The  conditions  for  stability  are  given  by  the 
Jury  test 


A  ( 1 )  >0 

m 

V1)  >  0 

and  1  >  |a_, 

j  m2 

Therefore,  the  parameters  must  satisfy 


aml  + 

am2 

aml 

am2 

and 

!  am2  j 

<  1 

(II . 14) 


J 

J 
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When  equations  (11.13)  and  (11.14)  are  satisfied,  the  impulse 
response  hm(n)  is  easily  shown  to  be  [10] 


h  (n)  = 
m 

z_1 

{ H 

m 

(z)  } 

a 

n 

sin (n+1 ) w 

m 

m0 

sin  co 

m 

Upon  substituting 

for 

z 

m 

★ 

and  z^  in  equation  (11.12)  and  equating 

A  (z  )  and  A  (z  ) ,  it  follows  easily  that 
m  m  m  m  1 


a  .  =  2r  cosu> 
ml  m  m 


a  _  =  -r 
m2  m 


(11.15) 


It  has  already  been  shown  that,  in  the  absence  of  noise, 
the  m'th  resonator  has  2  complex  conjugate  poles  located  within 
the  unit  circle  and  2  zeros  located  at  the  origin  of  the  unit 
circle  (equation  II. 11).  The  effect  of  the  addition  of  the  white 
noise  is  to  move  the  zeros  from  the  origin  towards  the  poles. 

When  the  noise  dominates  the  signal,  the  zeros  cancel  the  poles 
resulting  in  a  flat  spectrum. 

If  we  let  y  (n)  be  the  output  of  the  m'th  noisy  resonator,  then 

V">  -  +  T1 


(11.16) 


=  sm(n)  +  w' (n) 


where  w' (n)  is  white  with  variance 

2 


a 


w  1 


"hen,  the  power  spectrum  of  y  (n)  is 

m 
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Sy  (2)  =  °u  Hm(z)  2  +  <V 
1  m 


2  2 

amO°u  ,  _2 

- 1  +  °w ' 

A  ( z )  A  ( z  1 ) 
m  m 


amO°u  +  gw»Am(z)Am(z~1) 
Am(z)Am(z_1) 


If  we  let 


a2  B  ( z ) B  (z_1)  =  a2  n2  +  o2,A  (z)A  (z_1) 
nm'mv  mOu  w '  m  m 


(11.17) 


we  see  that  y  (n)  may  be  modeled  as  the  output  of  a  rational 

transfer  function  with  B  (z)  as  the  numerator  and  A  (z)  as  the 

m  m 

denominator  polynomials  with  n(n)  as  a  white  input.  Then 
- 1  -2 

B  (z)  =  1  -  b  ,z  -b0z  and  n(n)  is  a  zero  mean  white  sequence 
m  ml  mz 

2 

with  variance  o^.  It  is  known  [11]  that  a  Bm(z)  exists  which  has 
its  roots  on  or  within  the  unit  circle.  Equation  (11.17)  can  be 
written  as 

°n(1  -  bmlz~ 1  ~  bm2z_2)(1  ‘  bml2  “  bm2z2)  = 


(11.18) 


am.Oau  +  V  (1  -  amlz~ 1  "  am2z"2)  11  “  amlz  "  am2z2) 

Upon  expanding  the  above  equation  and  equating  terms  of  equal 
degree  on  both  sides  we  get 


22  2  , ,  2  2,  2,.  2  ,2, 
araO°u  +  V  (1  +  aml  +  am2  *  %(1  +  bml  +  bm2> 

aml°w*  (am2  '  U  -  bml%,bm2  '  1( 

2  _  2. 

m2  w  H  m2 


(11.19)  (b) 


Equations  (11.19)  (b)  and  (c)  imply 


aml(am2  "  1)bm2 


ml  am2 (bm2  '  1} 


(11.20) 


Equations  (11.19)  (a)  and  (c)  imply 


2  2  2  2  2 

a  „oz  +  o  ,(l  +  a  ,  +  a%) 
mO  u  w  ml  m2 


a  _,o  , 
m2  w 


1  +  b2  +  t/- 
ml  m2 


(11.21) 


m2 


or, 


k  +  H  +  am2  +  am2> 


m2 


2  2 

1  +  b,  +  bz_ 

ml  m2 


(11.22) 


'm2 


where 


X  = 


amO°u 


w 


If  the  noise  variance  a  ,  is  small,  then  X  >>  1.  Equation  (11.22) 


am2 


imDlies  therefore  that  - —  »  1  or  a  .  »  b  .  in  which  case 

bm2  m2  m2 


equation  (11.20)  shows  that  b^^  0,  Then,  for  hiqh  signal- 

to-noise  ratios,  the  zeros  are  located  near  the  origin  and  their 
presence  does  not  affect  the  spectrum  seriously.  The  exact 
locations  of  the  zeros  is  found  by  solving  equations  (11.20)  and 
(11.19)  simultaneously. 

II. 5  Perturbations  of  Pole  Positions 


The  addition  of  white  noiss  to  a  process  given  by  a  resonator 
mocel  tends  to  move  the  power  spectral  density  peak  and  to  broaden 
the  bandwidth.  It  has  been  seen  in  the  previous  section  that  this 
is  caused  by  the  zeros  being  shifted  from  the  origin  towards  the 
poles.  If  we  now  model  the  observation  again  by  an  all-pole  second- 
order  resonator  model,  the  positions  of  the  poles  will  have  changed. 


is 


It  is  this  perturbation  which  causes  an  expansion  in  the 
estimated  bandwidth.  In  the  following  this  effect  is  examined. 
The  output  of  the  m'th  resonator  is 


Ym(n)  =  s  (n)  +  w' (n) 


where  s  (n)  is  a  second-order  process. 


s  (n)  =  a  ,s  (n-1)  +  a  s  (n-2)  +  am.u(n) 
m  ml  m  m2  m  mO 


y  (n)  may  also  be  modeled  as  a  second-order  process, 
ym(n)  =  amlym(n-1)  +  am2ym(n-2)  +  amOu(n)‘ 

i  i  1 

In  the  absence  of  noise,  a  . ,  a  ~  and  a  .  are  identical  to 

ml  m2  mO 

aml '  am2  anci  am0  resPective^y -  The  autocorrelation  function 
r  (n)  of  the  process  y  (n)  has  been  shown  to  be  given  by  [10], 


ry  (n)  =  amOau  p  m  acos(lnK  “  V  +  °wl6(n) 


where , 


1  +  (- 


1  -  P 


1  +  P‘ 


1  +  o 


m,  2  .2 

■^r)  cot  ui 
2  m 


2  2  4 

l-o  1  -  2o  cos2ui  +  p 
m  m  m  m 


(11.23) 


=  arctan 


- r  cotu  . 

,  ,  2  m 

1  +  r 

m 


[II. 24) 


and,  p  and  o  are  related  to  the  parameters  a.  ,  a  by 
mm  ml  m2 


P  =  v-a  n 
m  m2 


CO  = 

m 


2 /-a 


arccos 
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1  a'  satisfy  the  Yule-Walker  equations  resulting  in 

aml ’  am2  "atl 


r  (0)r  (2)  -  r  (1> 

-  -P*« 

r2  (0)  -  r  U) 
ym  Ym 


(11.25) 


,  as  the  power  in  the  process 
Define  the  s ignal- to-noise  ratio  t  as 

s  (n)  divided  by  the  noise  power,  or 


r  (0)  -  ow, 

.. 


r  (0) 
3  m 


(11.26) 


=  Y  ~  1 

.  r  (0v  r  (1)  and  r  (2)  in  equation  (II. 25) 
Substitute  for  r  (0) ,  r  ym 

w  1 


111  .  1 

It  can  be  shown  that  (Appendix  A) 


,r  +  ^  +  c0 

r2  +  d. £  +  dn 


(11.27) 


where , 


cos2Wm  -  V  -  cos+mcosl2“i»  ' 

2cos2(um  -  V  -  3cos*mcos<2u,m  ‘  V 

cos2(%  -  >m  ' 

cos2!",  -  +m)  -  =°s*mcosl2“m  "  ♦*’ 

.  2(cos2*m  +  c°s2<»„  -  V  -  ®VOSP“-  ‘ 

.  3oos2»m  -  cos*n=°s(2%  -  ♦„)  +  COS2  (wm  -  *m> 
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The  distance 

pm 

from 

the  origin 

of 

the 

unit 

circle  to 

the 

location  of 

the 

pole 

is  related 

to 

the 

3  dB 

bandwidth 

B  in 
m 

Hz  of  the  m' 

th 

resonator  by  [  9  ] 

then , 


P 

P 


m 

m 


e 


e 


2 


B 

m 

2 


where  B  is  the  3  dB  bandwidth  in  Hz  of  the  noisy  resonator, 
ra  1 

Define  the  bandwidth  expansion  factor  (BEF)  as 


BEF  = 


m 


m 


=  1  + 


£n(% 

pm 


inp 


m 


(11.28) 


Equation  (11.27)  shows  that  the  BEF  depends  on  the  signal- to- 
noise  ratio  (SNR)  ,  £ ,  pm  and  wjn.  Families  of  curves  of  the  BEF 
as  a  function  of  the  SNR  in  the  range  -30  dB  to  +30  dB  have  been 
generated  for  3  different  frequencies,  .125  Hz,  .166  HZ  and  .25  Hz. 
They  are  shown  in  figures  (II. 4),  (II. 5)  and  (II. 6)  respectively. 

It  is  seen  from  the  graphs  that,  for  the  same  SNR,  the  BEF  is  an 
increasing  function  of  p  and  a  decreasing  function  of  u  .  For 
small  p's,  the  resonator  has  a  large  bandwidth  and  adding  white 
noise  (white  noise  has  flat  spectrum  between  -.5  Hz  and  +.5  Hz) 
does  not  affect  the  spectrum  as  much  as  for  large  p's. 

The  above  discussion  demons  crates  clearly  that  the  asymptotic 


resolution  of  an  all-pole  model  depends  in  a  non-linear  fashion 
on  both  the  bandwidth  and  frequency  of  the  noiseless  sianal  . 


Therefore,  in  general  an  increase  in  the  model  order  does  not 
uniformly  improve  the  spectral  estimates  of  all  the  spectral 
peaks  and  the  improvement  is  rather  unpredictable.  Curves  such 
as  given  in  figures  (II. 4)  -  (II. 6)  may,  however,  be  used  as 
a  general  guideline  for  the  amount  of  degradation  in  a  given 
application. 
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III.  ORDER  DETERMINATION  FOR  AR  MODELS 

The  previous  chapter  demonstrated  the  need  for  increased 
order  of  an  AR  spectral  estimator  to  improve  resolution  in  the 
presence  of  noise,  resulting  in  a  larger  variance  for  the 
estimates.  In  fact,  the  trade-off  between  the  bias  (or 
resolution)  and  variance  of  a  spectral  estimator  is  the  central 
issue  in  spectral  estimation  by  any  method.  For  the  traditional 
(Blackman  and  Tukey  type)  spectral  estimators,  this  trade-off 
is  reflected  in  the  choice  of  the  spectral  window  type  and  the 
maximum  lag  of  autocorrelation  function  used.  This  subject, 
referred  to  as  window  carpentering,  is  discussed  in  detail  by 
Jenkins  and  Watts  [12],  and  is  straightforward  because  resolution 
is  well-defined  in  terms  of  the  spectral  window  bandwidth. 

With  the  popularity  of  data  adaptive  (notably  the  auto¬ 
regressive  (AR) )  spectral  estimation  methods,  similar  resolution- 
variance  trade-offs  are  in  order.  Specifically,  well-defined 
methods  are  needed  to  determine  the  order  of  the  (AR)  spectral 
estimator  for  a  given  data  sample.  Furthermore  for  practical 
applications,  these  methods  need  to  be  on-line  and  as  much  as 
possible  objective  in  nature.  This  problem  is  complicated,  however, 
due  to  the  data  dependent  nature  of  the  resolution  of  the  AR  spectral 
estimator  as  shown  before  (e.g.,  no  well-defined  window  bandwidth). 
Therefore,  the  question  of  order  determination  for  the  spectral 
estimator  seems  to  be  best  posed  as  a  procedure  for  obtaining  a 
compromise  between  the  AR  model  fit  and  the  variance  of  the 
estimated  AR  parameters  as  a  function  of  the  model  order. 

Akaike  [13,14]  and  Parzen  [15]  have  recently  introduced  some 


methods  for  automatic  determination  of  orders  of  autoregressive 


processes.  One  method,  based  on  Akaike's  Information  Criterion 
(AIC) ,  has  gained  special  popularity.  In  this  report,  we  follow 
the  derivations  on  which  AIC  is  based,  introduce  appropriate 
modi fications  to  account  for  practical  estimation  procedures  and 
derive  a  new  information  criterion  designated  the  Conditional 
AIC  ( C A I C ) .  We  then  present  the  results  of  a  number  of  numerical 
simulations  that  compare  the  performances  of  AIC  and  CAIC  for 
moder  order  identification  and  spectral  estimation. 

III.l  Akaike's  Information  Criterion 

Akaike  derived  his  information  criterion,  AIC,  as  an 
estimate  of  the  asymptotic  relative  goodness  of  fit  of  the  model 
to  the  observation.  Although  his  derivations  were  based  on 
information  theoretic  arguments,  the  resulting  parameters  were 
the  same  as  the  maximum  likelihood  estimates.  In  this  section 
we  review  the  steps  involved  in  obtaining  AIC  as  they  pertain 
to  the  derivation  of  the  new  criterion.  We  assume  the  time  series 
to  be  described  by 
L 

x  =  Z  a.x.  •  +  u. ,  t=0,...N  (III.l) 

c  i=l  1  ^  c 


where  u  is  zero-mean  white  and  Gaussian  and  L  is  to  be  determined. 
Through  asymptotic  arguments,  Akaike  defines  an  information 
criterion,  related  to  the  maximum  likelihood  of  the  estimates  of 
ai  ,  a ^ ,  as: 

AIC (A)  =  (-2) In  (maximum  likelihood) 

+  E»  N  •  (  |  A  -  A  [  j  2 


(III .  2) 
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where  E^  denotes  asymptotic  expectation,  A  and  A  are  L  x  1  vectors 
of  the  coefficients  a^  and  their  estimates  a^.  The  practical 
AIC  which  is  related  to  the  full-information  likelihood  function 
of  a  Gaussian  process  is  then  given  by 

AIC(L)  =  N  In  (MLE  of  innovation  variance)  +  2L  (III. 3) 
and  the  order  L  is  chosen  that  minimizes  AIC(L). 


III. 2  The  New  Criterion 

Since  the  exact  maximum  likelih. ad  (full  information 
maximum- likelihood)  estimates  are  generally  not  available,  the 
conditional  MLE,  one  based  on  Yule-Walker  equations  or  Burg's 
algorithm,  of  the  innovation  variance  are  normally  used  in  (III. 3). 
We  propose  using  the  conditional  maximum  likelihood  (CML)  function 
in  (III. 2).  This  function  is  based  exactly  on  the  available  data 
and  we  believe  is  a  more  sensitive  indicator  of  the  behavior  of 
the  estimates  used  in  practice.  Thus,  in  the  following,  the  CML 
estimate  of  A  and  its  covariance  function  are  considered,  in  order 
to  obtain  tractable  expressions  for  (III. 2). 

The  conditional  (partial  information)  likelihood  function 
for  the  time  series  in  (III.l)  is  given  by: 


L  ( A ,  o  ^  ,x^, . . ,xL) 


N-L 


exp 


T 

-CDC 


( 2  7T  o  *■) 

u 


V] 


(in. 4) 


whe 


2  . 


re  au  is  the  variance  of  the  innovation  sequence  ufc, 


C1  =  [1  ,-aj,  ,-&2  >  •  •  •  »~aLJ  and 
^  =  =  k=L+l  Xk'i  +  1  Xk^+1  ‘ 


(III. 5) 
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Furthermore,  the  CML  estimate  of  oy  is  given  by: 


~  2  ~ 

a  =  C  DC/(N-L) 
u  -  — 


(III. 6) 


and  a  lower  bound  for  the  variance  of  the  estimates  of  a^  follows 
from  the  Fisher's  information  matrix  to  be 


var  [a,  ]  >  ^ 


(III. 7) 


where  A11  is  the  diagonal  element  of  the  inverse  of  the  (L  +  !)■ 


sample  covariance  matrix  of  xfc.  It  is  now  shown  that 


a  V1  -  1 

u  - 


(III. 8) 


The  innovations  variance,  ,  is  given  by  the  Yule-Walker 
relation  as 


Ju  ro  ~  [rl  •  -  *  rL^  -I,1 


(III. 9) 


lrLJ 


3 


where  AL+1  is  the  L+l-st  order  covariance  matrix  given  by: 

....  rT 


-  L+l 


L-l 


r  ....  r 
L  o 


The  first  diagonal  element  of  the  inverse  of  AT  -  ,  is  then 

given  by : 


ill 


'-L+l  ■ 


(III . 10) 


Formula  (III. 9)  can  be  rewritten  as 


f 
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v 


Resulting  in 


a 


2 


u 


(III. 13) 


as  required. 

We  now  proceed  to  define  an  expression  for  (III.  2)  based  on 

2 

the  CML  estimates  of  A  and  au  .  The  expression  for  the  CML  given 
in  (III. 4)  is  now  substituted  in  (III. 2)  for  the  maximum  likelihood 
and  using  (III. 8)  for  the  second  term  in  (III. 2)  and  (N  -  L)  for 
N  we  have : 


CAIC(L)  =  (N  -  L) ln(2nou2)  +  (a  -  1)L  (III. 14) 

The  factor  a  >  1  is  included  to  account  for  the  asymptotic  nature 
of  the  criterion  and  the  fact  that  (III. 8)  is  a  lower  bound  for 
the  variance  of  a,.  A  similar  parameter  was  also  suggested  for 
AIC  [16]  and  in  [14]  Akaike  discusses  a  possible  approach  for 
choosing  a.  Since  CAIC(L)  as  given  by  (III. 4)  is  dependent  on 
the  variance  of  xfc,  the  test  is  standardized  by  introducing  a 
normalized  innovation  variance  so  that 


CAIC(L)  =  (N  -  L) In [au2/ (var  xfc) ] 
+  (a  -  1)  L 


(III. 15) 


Thus  CAIC(O)  =  0.  The  factor  a  is  chosen  to  give  more  or  less 
weight  to  the  error  in  the  estimation  of  the  parameters.  In 
other  words,  resolution  can  be  increased  at  the  expense  of  the 
variance  of  the  estimates  by  decreasing  a.  We  have  found, 
empirically,  values  of  3.5-4  to  give  the  most  stable  and  reasonable 
indication  of  the  order. 


31 


III. 3  Simulation  Results 

We  have  tested  the  performance  of  CAIC  relative  to  AIC  on 
a  number  of  time  series  models  reported  previously.  The  data 
included  normal  as  well  as  uniform  distributions.  The  estimates 
were  based  on  CML  (least-square)  and  Yule-Walker  methods.  In 
the  great  majority  of  cases,  CAIC  performed  as  well  or  superior 
to  AIC.  Examples  of  these  can  be  found  in  [17] .  Some  estimated 
spectra  based  on  orders  determined  by  AIC  and  CAIC  are  also  shown 
in  figures  (III.l)  -  (III. 3).  Yule-Walker  equations  with  auto¬ 
correlation  function  estimates  given  by 


i  N:£ 

ri  N  XjXj  +  i 


were  used.  The  example  shown  in  Figure  (III.l)  indicates  the 
relative  stability  of  CAIC.  Figure  (III. 2)  shows  that  the  model 
order  chosen  by  AIC  results  in  spurious  peaks,  while  giving 
higher  peak  resolution  than  the  CAIC  based  one.  Figure  (III. 3) 
shows  that  an  increase  in  white  noise  level  increased  the  AIC 
order  to  the  point  that  spurious  spectral  peaks  became  pronounced 
while  CAIC  remained  nearly  the  same,  showing  the  relative  stability 
of  CAIC. 
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FIGURE  III.  3  Log  power  spectrum  of  xt=sin  ( .  6iTt+30 °  )  +0 . 707sin  ( .  27it+60 °  ) 
+nt;  an2=-316'  N=100;  CAIC=5 ,  AIC=15. 
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IV.  SOME  SUBOPTIMUM  ARMA  SPECTRAL  ESTIMATORS 


As  shown  in  Chapter  II,  the  presence  of  additive  noise  on 
the  observations  of  an  AR  process  implies  that  they  are  an 
ARMA  process.  A  popular  approach  to  circumvent  the  degradation, 
wrought  by  this  model  change,  of  the  AR  spectral  estimate  is 
to  derive  spectral  estimators  based  on  the  ARMA  model.  Box 
and  Jenkins  (  7  ]  among  others,  provide  an  algorithm  to  give  a 
close  approximation  to  the  maximum  likelihood  (ML)  estimate  of 
the  ARMA  parameters,  assuming  Gaussian  observations.  This  is 
considered  the  optimum  estimator  in  view  of  the  desirable 
properties  of  maximum  likelihood  [18].  But  optimum  estimators 
have  computational  disadvantages,  described  in  the  next  chapter, 
which  have  compelled  many  researchers  to  consider  alternative 
criteria  that  yield  spectral  estimators  with  greater  computational 
efficiency.  There  is  presently  considerable  activity  in  this 
area,  appropriately  termed  suboptimum  ARMA  spectral  estimation. 

In  this  chapter  a  new  suboptimum  scheme  is  reported  for 
estimating  the  power  spectral  density  (PSD)  of  an  ARMA  process 
of  known  orders.  After  a  preliminary  data  reduction,  this  scheme, 
called  the  least  squares  (LS)  estimator,  minimizes  a  sum  of  squared 
quadratic  functions  of  the  AR  coefficients  using  a  nonlinear  least 
squares  algorithm.  The  poles  of  the  estimated  PSD  are  found  from 
the  minimizing  AR  coefficients,  and  zeros  are  found  from  quadratic 
functions  of  these  coefficients.  Note  that  these  results  have 
been  published  [ 19] . 

The  general  idea  of  least  squares  fitting  the  ARMA  parameters 
is  not  new,  and  various  other  approaches  have  been  suggested 


i  -sV  l  ;  ivt<4 
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(see,  e.g.,  [7]  and  [20]).  The  scheme  discussed  here  is. 


however,  analogous  to  a  minimum  mean-squared  error  estimation 


of  the  parameters  appearing  in  the  estimator  discussed  in  [21] 


and  [22].  A  modification  of  the  latter  estimator  that  is  based 


on  the  modified  Yule-Walker  equations  (the  MYW  estimator)  is 


also  presented,  in  which  the  problem  of  negative  excursions  of 


the  estimated  PSD  is  corrected  by  tapering  the  estimated  moving- 


average  autocorrelation  function.  Examples  are  shown  that  compare 


the  performance  of  these  ad  hoc  techniques  to  the  approximate 


maximum  likelihood  method  of  3ox  and  Jenkins. 


IV.  1  The  Spectrum  of  an  ARI1A  Process 


Assume  that  we  observe  xt,t=l,...,N  where  xfc  is  stationary 
and  Gaussian  of  mean  zero,  and  that  xfc  fits  an  ARMA  (L,M)  model. 


Then  we  can  write 


aixt-l  =  ut 


(IV. 1) 


where  ai  are  the  AR  parameters  for  the  ARMA  model  and  ufc  is  the 


MA  residual  sequence  given  by 


ut  =  et  “  l  bict-i 


(IV. 2) 


where  are  the  MA  parameters  and  is  a  zero-mean  uncorrelated 


normal  sequence  of  variance  Define 


A  [1,— a^,«..,— a^]. 


Then  the  PSD  of  x^_  is  given  by 


-  SuU)  '2 


(IV. 3) 


T  —  1 

where  =  [l,z  ,...,z  ],  z  being  the  z-transform  operator 
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1 


c  f  x  is  the  PSD  of  ut.  Writinq 
z  ,  e3'^,  and  where  Su(z)  is  tne  t 

-M  .  ,Mv.  (IV. 4) 

2„  +  o  (.-1  +  z)  +  •••  +  6M(Z  +  Z  )] 

Su(z)  =  ou[l  + 

we  can  egress  -  .SO  of  ut  in  ter„s  of  if  •>* 

M  noticed  -  ""  COefflC7  0 pp. 

of  ot  ,1~  ^  -  conditional  elation 

218-225] 


E(utiut.i)  -  3iut-i 

,  3  is  the  best  linear  predictor  of  ut  aiven 

and,  in  general,  %  is  tne 

mpan-sauare  error  (mmse)  sense, 
in  the  minimum  mean  squai 

t-i 

-o  mhe  LS  Estimator 

— -  4_„  ft  is  obtained  by  mminu^  . 

A  least  squares  estimate  o  t 

given  by 
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(IV. 5] 
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t=L+i+l 


(Alt,t-0)  -  h^ht-i.t-L-ii 


where 


li,  jl 
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The  derivative  with  respect  to  6,  vanishes  for 


ATRiA 


(IV. 6) 


where 
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N  T 

Ri  =  t=L+i+l  'Xtt-t-L)X  (t-i.t-L-il 


+<XIt,t-L|xT[t-i,t-L-i]'Tl 


R°i  t=L+i+l  [XIt-i,t-L-i]X  [t-i,t-L-i] 1  * 


That  0^  minimizes  is  seen  by  taking  the  second  derivative 


2ATR  .  A 
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[t-i ,  t-i-L] 


>  0. 


(IV. 7) 


This  also  establishes  a  nonzero  denominator  for  0^  in  (IV. 6). 

To  get  an  estimate  of  the  variance  of  ufc  we  use  the  sample 
moment 


(N-L)u  = 


N 

£ 

t=L+l 


u,  = 


N 

1 

t=L+l 


A^X 

A  X[t,t-L]X 


[  t ,  t-L] ' 


T 

A  R  A 
oo 


(IV. 8) 


where  R  is  given  in 
The  fact  that  ufc 
one  method  to  fit  the 
equations  (3^  =  0,  i  = 
that  the  denominators 
simpler  system 

ATRiA  =  0, 


(IV. 6)  . 

is  MA(M)  implies  that  ( 
vector  A  to  the  data  is 
M  +  1 ,  ...  , M  +  L  from 
are  nonzero,  we  can  ins 

i  =  M  +  1,  ...  ,M  +  L. 


^  =  0  for  i  >  M,  so 
to  solve  the  set  of 
(IV. 6).  Recognizinq 
tead  solve  the 

(IV. 9) 


The  locus  of  solutions  for  each  of  the  above  equations  is  a  quadric 
surface  (24,  pp.  287-294],  and  there  is  no  solution  to  the  system 
if  the  surfaces  do  not  jointly  intersect  at  least  one  point.  In 
order  to  obtain  a  value  for  A  whether  or  not  there  is  a  common 
intersection,  we  choose  A  to  minimize  Q  given  by 
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Q  =  E  (A  R .  A)  .  .  (IV. 10) 

i=M+l  1 

The  value  of  Q  is  zero  at  any  point  of  intersection  for  all  the 
surfaces  in  (IV. 9).  The  vector  A  obtained  via  minimization  of  Q 
will  henceforth  be  referred  to  as  the  LS  estimate.  Using  the 
LS  estimate,  we  next  find  8^,  i=l,  ...  ,M  using  (IV. 6)  and  ou 
using  (IV. 8)  and  substitute  these  values  in  (IV. 4)  and  (IV. 3)  to 
get  an  estimate  of  the  PSD  of  xfc,  henceforth  termed  the  LS 
spectral  estimate. 

IV. 3  The  MYW  Spectral  Estimator 

The  rational  spectral  estimator  in  [  j  (the  MYW  estimator) 
produces  A  to  solve  the  modified  Yule-Walker  equations,  which 
for  an  ARMA  (L,M)  process  are 


39 


and  the  estimated  PSD  is 


Sx(z)  =  [ru(0)  +  ru(l) (z-1+z) 


+ . 


r  (M)  (z 
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(IV. 14) 


Note  that  the  LS  and  MYW  estimates  of  A  are  asymptotically 
equivalent.  To  see  this,  observe  that 


P 


lim (N- L-k) 
N-’-00 


+ 


A 


T 

k 


(denotes  convergence  in  probability)  [25]  where  A^  incorporates 
the  mth  lag  autocorrelation  rx(m)  of  xfc  according  to 


13 


rx(k+i-3)  » 


i  —  1,  ...  ,  L+l ,  j  ■  1,  ...  ,  L+l  , 

and  that  p  lim^^C^  =  ||Aij||  ,  i  =  l,  ...  ,L;  j  =  l,  ...  ,  L+l.  The 

k 

solution  to  ( IV .11)  obtained  by  replacing  CM+^  by  its  limiting  value 
also  solves  (IV. 9)  with  (N-L-k)  replaced  bv  its  limiting  value. 

There  is  one  other  solution  to  this  asymptotic  form  of  (IV. 9), 
but  we  conjecture  that  it  is  outside  the  stationary  region  for  A. 

The  spectral  estimate  given  by  (IV. 14)  is  not  guaranteed  to 
be  a  nonnegative  function  of  z  =  e-^71^  for  fe  [0,1/2].  For  instance, 
if  observations  consist  of  two  additive  narrow-band  AR(2)  signals 
having  center  frequencies  in  close  proximity,  the  true  z-spectrum 
consists  of  two  closely  spaced  poles  just  inside  the  unit  circle 
and  a  zero  just  inside  the  unit  circle  and  between  the  poles  in 
frequency.  An  error  in  estimating  the  numerator  polynomial  in 
(IV. 14)  can  cause  what  should  be  a  near-zero  positive  value  at 
the  bottom  of  the  trough  in  the  frequency  response  of  the  numerator 
to  be  a  near-zero  negative  value,  thus  making  the  PSD  estimate 
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negative  in  a  region  near  the  peak.  An  effective  counter  is 
to  reduce  the  depth  of  the  trough  by  multiplying  the  r  ( i)  by  a 
taper  which  slightly  reduces  the  frequency  resolution  of  the 
estimated  numerator.  Use  of  the  linear  taper 

Ti  =  1  -  i/K,  i=0 , 1 ,  ...  , M ,  K  >  M  (IV. 15) 

has  been  successful  in  eliminating  negative  excursions  of  the 
PSD  estimate  when  used  to  produce  ru(i)t  according  to 

vi}t  =  Vu(i)'  i=0'  •••  'M- 

Another  strategy,  used  in  [20],  is  to  replace  the  numerator  in 
(IV. 14)  by  the  periodogram  of  ufc.  For  short  data  records,  however, 
the  periodogram  often  cannot  adequately  represent  the  movinq 
average  spectrum,  resulting  in  inaccurate  indication  of  the  power 
under  the  peaks  of  the  ARMA  spectrum.  In  tests  of  the  MYW 
estimator,  all  negative  excursions  were  eliminated  for  K  >>  M. 

No  negative  excursions  were  noted  for  the  LS  algorithm,  but  the 
existence  of  a  guarantee  has  not  been  investigated.  It  is 
emphasized  that  the  use  of  the  numerator  taper  does  not  sacrifice 
the  resolution  of  the  ARMA  spectral  estimates.  This  is  due  to 
the  fact  that  resolution  is  mainly  determined  by  the  poles  and 
the  MA  spectra  of  interest  are  relatively  smooth. 

IV. 4  Simulation  Results 

Figure  (IV. 1)  shows  LS  estimates  of  the  PSD  of  20  realizations 
of  xfc  =  wfc  +  Yfc  +  0.5  nfc  where  wfc  =  0.4  wfc_^  -  0.93  wt_2  +  and 
yt  =  -0.5  “  0.93  yt_2  +  nt  and  where  nt»  and  nfc  are 

mutually  independent  i.i.d.  Gaussian  of  mean  zero  and  unit  variance. 
Q  in  (IV. 10)  is  minimized  using  the  conjugate  gradient  technique  [26] 


with  A  starting  at  the  origin.  The  combined  effect  of  the  plots 
is  to  suggest  a  bias  in  the  estimator  which  smears  the  spectra. 

In  Fig.  (IV. 2)  the  tapered  MYW  estimate  is  depicted  for  the  same 
set  of  realizations.  K  in  (IV. 15)  is  set  to  the  minimum  value 
that  succeeds  in  eliminating  all  negative  excursions  of  the 
spectral  estimate.  In  this  set,  the  untapered  MYW  estimates  of 
seven  realizations  went  negative,  and  K  =  34.  Figure  (IV. 3)  shows 
the  uncondi tional  least  squares  Box  and  Jenkins  estimate  for  the 
same  set  of  realizations.  It  is  apparent  that  the  performance 
of  the  MYW  estimator  is  nearly  as  good  as  that  of  the  Box  and 
Jenkins  estimator  for  most  realizations.  The  disturbing  tendency 
of  the  MYW  estimator  to  produce  an  occasional  aberrant  estimate 
is  also  exhibited  in  the  figure. 

The  simulations  provide  evidence  that  the  computational 
simplification,  compared  to  the  optimum  schemes,  provided  by 
data  reduction  in  the  LS  estimator  is  accomplished  at  the  expense 
of  a  substantial  tradeoff  in  statistical  efficiency.  The  results 
also  suggest  that  the  MYW  estimator  has  higher  statistical  efficiency 
than  the  LS  estimates.  These  observations,  coupled  with  the  high 
computational  efficiency  of  the  former,  lead  us  to  conclude  that 
the  MYW  estimator  is  the  superior  of  the  two  suboptimum  spectral  . 
estimators  considered. 
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Fig.  IV. 1.  LS  spectral  estimate,  N  =  125. 
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Fig.  IV. 3.  Unconditional  least  squares  Box  and  Jenkins 
spectral  estimate,  N  =  125. 
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V.  STATISTICAL  CLASSIFICATION  OF  SOME  ARMA  SPECTRAL  ESTIMATORS 


Many  of  the  spectral  estimates  appearing  in  the  engineering 
and  statistics  literature  over  the  last  two  decades  fall  into 
one  of  two  classes.  In  the  minimal  sufficient  (MS)  class,  the 
AR  and  MA  parameters  are  adjusted  simultaneously  to  give  a  least 
squares  or  approximately  least  squares  fit  to  the  observed  data 
vector.  The  name  "minimal  sufficient"  acknowledges  that  no 
reduction  of  the  data  is  a  sufficient  statistic  for  stationary 
Gaussian  time  series  having  zeros  in  the  power  spectral  density 
(PSD)  [31].  The  resulting  estimates  are  approximately  maximum 
likelihood  (ML)  for  Gaussian  data.  In  the  sequel,  it  is  assumed 
that  the  data  is  a  Gaussian  time  series,  and  the  ML  ARMA 
parameter  estimate  is  referred  to  as  optimum,  for  reasons 
discussed  in  Chapter  IV.  Thus,  optimum  solutions  are  contained 
in  the  MS  class.  The  equations  for  the  least  squares  solution 
are  highly  nonlinear,  and  thus  estimators  in  this  class  are 
typified  by  nonlinear  optimization  algorithms  and  other  iterative 
approaches  and  their  concomitant  pitfalls,  i.e.,  nonconvergence 
or  convergence  to  local  rather  than  global  extrema,  the  need  for 
preliminary  ARMA  parameter  estimates  to  start  the  iterations, 
and  large  computer  memory  and  time  requirements.  The  auto¬ 
correlation  function  ~iass  is  composed  of  estimators  which 

utilize  a  fixed,  finite  number  of  lags  of  the  sample  ACF .  The 
estimation  equations  are  usually  linear  or  quadratic,  and  often 
the  AR  parameters  are  estimated  alone,  after  which  estimates 
either  of  the  MA  parameters  or  of  some  function  related  to  the 
PSD  of  the  MA  residuals  are  computed  using  the  AR  estimates.  All 


of  the  estimators  in  this  class  are  suboptimum,  the  truncated 
sample  ACF  being  a  data  reduction.  Further,  an  optimum  fitting 
requires  that  all  parameters  be  adjusted  simultaneously. 
Suboptimality,  then,  is  the  price  paid  for  the  gains  in 
computational  simplification  and  relaxed  memory  requirements 
of  estimators  in  this  class. 

In  its  emphasis  on  computational  efficiency,  the  literature 
has  left  largely  unanswered  questions  regarding  the  statistical 
efficiency  of  estimators  in  the  ACF  class.  We  provide  in  this 
chapter  an  evaluation  of  Fisher's  information  matrix  for  the 
truncated  sample  ACF  of  ARMA  processes.  Only  the  asymptotic 
(infinite  observation  record,  i.e.,  N  -*■  co)  case  is  treated  to 
make  the  analysis  tractable,  and  the  magnitude  and  angle  of  poles 
and  zeros  of  the  PSD  are  taken  as  parameters  of  interest.  The 
results  will  yield  asymptotic  bounds  on  the  statistical  efficiency 
of  any  estimator  in  the  ACF  class. 

V.l  Classification  of  Some  ARMA  Spectral  Estimators 

In  keeping  with  the  goal  of  obtaining  asymptotic  results, 
the  estimators  to  be  discussed  are  classified  according  to  the 
limiting  form  that  the  statistic  on  which  they  are  based  takes 
as  N  approaches  infinity.  The  MS  class  includes  the  methods  of 
Tretter  and  Steiglitz  [32],  Hannan  [33],  Akaike  [34],  Konvalinka 
and  Matousek  [35],  and  Box  and  Jenkins  [7,  pp.  231-235].  The 
methods  of  Walker  [36],  Hsia  and  Landgrebe  [37],  Graupe , 

Krause  and  Moore  [38],  Sakai  and  Arase  [39],  Satorius  and  Alexander 
[40],  Kaveh  [21],  Kinkel  et  al.  [22],  and  Bruzzone  and  Kaveh  [19] 
are  induced  in  the  ACF  class.  Cadzow's  method  [20]  is  based  on 
a  reduction  of  nearly  all  N  lags  of  the  sample  ACF  and  hence  does 
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not  fit  into  either  class.  We  have  observed  a  tendency  toward 
large  frequency  errors  in  Cadzow's  method  as  well  as  an 
exaggeration  of  the  sharpness  of  spectral  peaks  while  requiring 
nearly  as  much  computer  time  and  memory  as  methods  in  the  MS 
class.  These  facts,  in  conjunction  with  the  suboptimality  of 
the  method  as  a  result  of  the  use  of  a  data  reduction  as  well 
as  non-simultaneous  determination  of  the  AR  and  MA  parameters, 
lead  us  to  dismiss  Cadzow's  method  from  further  analysis.  A 
number  of  estimators,  e.g.,  [36],  which  are  based  on  the  sample 

ACF  yet  claim  to  be  asymptotically  efficient  from  a  subclass  of 
the  ACF  class.  There  is  no  discrepancy  here,  however,  in  that 
efficiency  obtains  in  the  limit  as  first  the  number  of  observations, 
then  the  number  of  sample  ACF  lags,  approach  infinity.  In 
practice  we  are  not  at  liberty  to  extend  the  estimators  in  such 
fashion,  so  these  estimators  are  analyzed  assuming  a  finite 
truncation  of  the  sample  ACF. 

Looking  more  closely  at  the  ACF  class,  we  see  that  [37], 

[21],  [22],  [40]  and  [19]  first  compute  an  estimate  of  the  AR 

parameters,  and  then  use  these  to  estimate  some  function  related 
to  the  MA  parameters.  Thus,  the  statistical  efficiency  of  these 
estimators  is  determined  largely  by  the  efficiency  of  the  AR 
parameter  estimates  in  the  first  stage.  The  AR  parameters 
determine  the  pole  locations  of  the  ARMA  PSD,  and  consequently 
they  have  far  greater  influence  on  its  shape  than  do  the  MA 
parameters  in  the  problem  of  high  resolution  spectra.  It  is 
clear  that  the  asymptotic  efficiency  with  respect  to  only  the 
estimated  pole  locations  of  various  estimators  is  a  meaningful 
criterion  by  which  to  judge  them.  For  example,  [37],  [21],  [22], 


[40]  and  any  method  based  on  modifying  the  numerator  portion 
of  their  estimated  PSD's  (as  in  Kay  [41] )  are  equivalent  by 
this  criterion,  a  result  that  agrees  with  our  experience. 


The  two  popular  sample  ACF's 


1  ^ 

Cx(l)  N  tf1  xtxt+i  ;  1 


and 


<V.l] 


1  N_i 

fPT  t^1  xtxt+i  '  1  0,1,‘'''k 


are  asymptotically  equivalent  for  finite  k.  We  take  as  the 

T 

statistic  of  interest  C^k  k  ^  =  [cx  (k^)  ,cx  (k^+1)  , . . .  ,cx (k2) ]  , 

k^  <  ^2"  T^i-S  allows  us  to  take  into  account  the  effect  of  not 
using  low-lag  ACF  values,  as  is  the  case  in  the  five  estimators 
listed  in  the  previous  paragraph,  wherein  cx(o)  is  not  used  to 
estimate  the  AH  parameters.  The  asymptotic  Fisher’s  information 
matrix  is  derived  for  any  subset  of  the  poles  and  zeros  by 
appropriate  choice  of  the  parameter  vector. 

We  mention  briefly  the  work  of  Gersch  [42],  who  provides 
the  asymptotic  covariance  matrix  of  the  AR  parameter  estimates 
gotten  from  the  modified  Yule-Walker  operations,  and  that  of 
Sakai  and  Tokumaru  [47],  who  give  the  covariance  of  the  estimated 
power  spectrum  gotten  from  using  [21]  or  [22]  (these  are 
equivalent) .  Little  else  has  been  done  to  analyze  the  suboptimum 
ARMA  spectral  estimators,  and  while  these  results  are  useful, 
they  lack  the  general  applicability  of  the  analysis  we  now  undertake 


rr 
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V.2  Evaluation  of  Fisher’s  Information  Matrix 

We  begin  by  evaluating  the  asymptotic  covariance  matrix  of 
the  sample  ACF  for  an  ARMA(L,M)  process,  where  we  assume  L  >  0, 
using  the  notation  and  assumptions  of  Chapter  IV.  It  is  known 
that  (12,  p.  181],  for  a  stationary  time  series  xfc , 

1  °° 

dm  cov  (c  (k)  ,c  ( l)  ]  =  rr  I  y  (i)  y  (i+t-k) +y  (i  +  2,)  Y  (i-k)  (V.2) 

j_5=—00 

rfhere  y(i)  =  E(xfcxt+j,).  The  ACF  of  an  ARMA(L,M)  process  is 
recursive  beginning  with  lag  M.  Starting  values  are  given  by 
the  appropriate  system  below  {see  (44]), 


1  -aL  ~a2  ~aL~l  "aL 

Y  (0)  ' 

r° 

0 

^"*a2  ~a3  *  '  *  ~aL  ® 

Y  (1) 

“a2  -aI_a3  a4  •  •  •  0  0 

•  *  *  *  * 

Y  (2) 

• 

=  o  2 
z 

D 

2 

- 

•  •  •  *  • 

y(L-l 

• 

~aL  ~aL-l  ~aL-2  '  '  •  'ai  1 

Y  (L) 

D 

W 

;  L  >  M  >  0 
(V. 3a) 
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0 

0 

0 

...  0 

Yd) 

D 

1 

~a2 

"al_a3 

L"a4 

...  0 

0 

0 

0 

...  0 

Y(2) 

D 

2 

~a  L 

_aL-l 
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(V. 3b) 
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where 


D  =  a  2  r  V  (a  -bt)B  detlo,  )  *  ot*Bj ,  (V.4) 

3=1+1  H=1  J  J 

with  =  1,  B.  =  for  =  0  for  i  >  M,  a^  =  0  for 

3 

i  >  L,  and  1=0  for  j  <  i.  Also, 
i 


ak  = 


al  a2 


- 1  a 

0  - 


1  **•  k-1 

1  •••  ak- 2 


0  0 


;  k  >  1 


(V.  5) 


;  k  =  0 


where  a.  =  0  for  i  >  L.  For  M  >  1  and  L  >  1  we  extend  the  vector 
of  starting  values  to  lag  M+L-l  using  the  AR  recursion 


Y (k)  =  1  a. y  (k-i) . 

i  =  l  1 


( V .  6 ) 


Although  this  recursion  can  be  used  to  extend  the  ACF  arbitrarily 
far,  the  infinite  sums  in  (V.2)  are  simplified  by  expressing  the 
ACF  in  terms  of  its  poles,  i.e.. 


Y(i)  =  ~  ^G.1  M  ,  i  >  M  , 

j  =  l  1  3 


( V .  7  ) 


where  the  are  zeros  of  the  characteristic  equation 


1  -  1  a .  z  =  0  , 

i=l  1 


(v.  8: 
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and  Aj  is  the  residue  of  G ^ .  Usinq  the  extended  ACF,  we  find 
the  Aj  by  solving 


1  ...  1 


y  (M-t-1) 


(V .  9 ) 


G , L_ 1  G0L-1  ...  GtL  1  XT  y  (M+L- 1 ) 

i  i  ^  ij  u 


The  infinite  sums  in  (V.2)  can  be  expressed  in  terms  of  one-sided 


infinite  sums,  as 


2  I  y2  (i)  +  y2  (0)  ;  j  =  0 


Z  Y(i)v(i  +  j)  = 


2  I  Y  ( i )  Y  ( i  +  j )  +  I  Y(i)Y(j“i)  ;  j  odd  (V.10) 

i=0  i=l  j 


2  Z  y  ( i ) Y (i  + j )  +  I  Y  ( i ) Y  ( j“i )  +  y  ;  j 


and  the  infinite  sums  can  be  separated  into  terms  of  the  ACF 
directly  influenced  bv  the  MA  parameters  and  those  followinci  the 
recursive  form  (V.7),  as 


Z  y  ( i  )  Y  ( i  +  j  )  =  l  y(i)Y(i  +  j)  +  Z  Y(i)Y(i  +  j) 
i=0  i=0  i=M 


But,  from  (V.7) , 


to  co  l  L  ... 

-  Y  ( i )  Y  ( i  +  j  )  =  Z  Y.  Z  \\ „Gm1"MG*  +  :,"M  = 


i=M  m=l  n=l 


m  n  m  n 


L  L  «  r  T  r  3 

z  ::  X  A  Z  _  i-M_  i+i-M  “  “  ,  ,  _ 

m=l  n=l  i=M  m  n  ,  .  m  n  1-G  G  J 

m=l  n=l  m  n 

(V.ll) 


51 


where  the  assumption  that  all  poles  are  inside  the  unit  circle 
allows  the  interchange  of  summation  and  guarantees  convergence 
to  the  closed  form  result. 

Substituting  (V.ll)  into  (V.10),  we  have 


,  L  L  XX 

y2(0)  +  2  I  I  Irf-g-  ;  1=0 

m=l  n=l  m  n 


CO 

I  y  (i)  y  (i+j ) 
i=-°° 


M-l 

j-1 

2 

E  y ( i ) y (i+j ) 

+  E 

i=0 

i=l 

L  L 

Gnj 

2  E  E  X  X 

t  ,  m  n 

l-G  G 

m=l  n=l 

m  ] 

M-l 

i-i 

E  y(i)y(i  +  j) 

+  E 

i=0 

i=l 

7  ■  L 

L 

Y  (J)  +  2  E 

m=l 

E  X 

n=l  n 

n 


( V .  12  ) 


j  even. 


m  n 


Finally,  we  note  that  (V.2)  can  be  rewritten 


£im  cov[cx(k)  ,c  U)  ]  =  ^  E  y  ( i )  y  ( i+2.-k)  +y  (i)  y  (i- l-k)  ,  (V.13) 

N-VOO  *  —  oo 


so  that 


|U2<°)  *  2  S M  *  2  art  n-l^n]' 


£  £  Y 

N  i  =  l 


2(i)  +  z  7  (i)-Y  (i+2?-)  + 


E  Y(iW(2l-il  + 
i  =  l 


L  L  ,  n  S  \  + 

+  r,  „L  V»'PS 

m=l  n=I  -1 


i.(v2(0)+x2(2))  ;  I  =  k  >  0 

N 


cov  [c  OO  >c%a)  1 


2  j^1  Y  (i)  (y  (i+*+k)+Y  + 


dim  covlcx 


(V.14) 


N  1  =  0 


d+k-l 
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2,-k-l 
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2  ...  y  V  (i)  Y  (8,-k-i)  + 
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L  L  ,  ,  n  n  _} 

r1  A  x»  " 

m=l  n-1 
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where  terms  used  in  the  finite  sums  are  obtained  using  (V.6). 
Note  that  the  above  assumes  i  >  k ,  £  >  0  and  k  >  0.  This 
presents  no  problem  in  evaluating  the  covariance  matrix,  for 
cx(k)  =  cx(-k)  implies  that  negative  values  of  l  and  k  are  not 
needed,  and  the  toeplitz  form  of  the  covariance  matrix  implies 
that  cov[c  (k)  ,c  (*.)  ]  =  cov  [c  ( l)  , c  (k )  ]  .  The  case  L  =  0  is 

A  X  XX 

handled  by  (V.3b)  with  (V.4)  and  (V.5).  (V.10)  is  modified  to 
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Then,  using  (V. 2) , 

we 
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i ,  for  L  =  0 , 

Y (i) Y ( j-i) 
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(V. 16) 


/  3 , 5  ,  .  .  . 


+ 


1 

Y  (i  )  Y  (  2-k-i ) 


;  2-k=2 ,4,6,... 
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The  probability  density  function  (p.d.f.)  of  the  sample  ACF 
of  an  ARMA  process  is  asymptotically  multivariate  Gaussian  [  ] 

If  we  denote 


F  [k  k  ]  =  Iy  (kl>  ' Y  (ki+l)  '  •  •  •  'T  (k2)  1  ;  k2  >  kl 


Atk1.K-2I  =  ““  ElC(k1.)-.2rr[k1,k2l),C[|<i.k2i'rIK1,k2],T  '  {V .17) 


^ 2  J  L  l '  2  J 


then  the  asymptotic  p.d.f.  of  Cr,  ,  ,  is  given  by 

l-Kl  'K2  J 


k2~kl 


lim  f  (Cr,  ,  ,  |  G ,  P ,  a  4 )  =  (2n  ] 

N-oo  lkl'k2J  e 


[k1<k2]  1  exp  "  2(C[klfk2] 


[k1,k2})  A[kL,k2]  lu[k1,k2]~1  [k1,k2] 


T.ri 


(V. 18) 


where  G  represents  the  vector  of  poles  and  P  the  vector  of  zeros. 
Fisher's  information  matrix  is  given  by  [45]  (denote  lim  f  =  f^) 

T 

vv  -  ~e  \(['Vn£(C|ki'i!2]|G'p'aE  n  1  ■  (V-19) 


We  take  ‘jQ  ^  I  G  ^  1  ,  •  -  .  ,  i^l/2  7  ^  '  ’  •  •  ' ^  l/ 2  f  ^1  9  *  *  *  r  ^*M/2  9  ^  1 9  *  *  *  1  ^M/  2  J  9 
where  the  assumption  that  poles  and  zeros  occur  in  complex  conjugate 
pairs  halves  the  size  of  the  parameter  vector.  The  ate  angles  of 
the  pole  pairs,  and  zeros  P^  solve 


1  -  b ^z 


,  -M  n 
••  ~  bMz  =  0  ' 


(V. 20) 


each  pair  having  angle 


I 
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We  have  been  unable  to  obtain  the  derivatives  in  (V.19) 

due  to  the  complicated  nature  of  the  p.d.f.,  and,  even  if  they 

were  available,  it  is  doubtful  that  the  expectation  would  yield 

to  a  concise  solution.  Numerical  methods  are  not  well  developed 

for  multiple  integrals,  and  these  have  the  added  complication 

of  being  improper,  so  we  consider  Monte  Carlo  methods.  The 

most  basic  approach  is  to  generate  several  thousand  points 

(Cr,  ,  ,) •  distributed  according  to  (V.18)  and  to  approximate 

1 ' K  2  1 

the  expectation  by  a  sample  mean 


W 


1=1  O  O  1  2 


1  |  h  2  ^ 

g  K,  U7e  If  ((cfk  .  ,).  G ,  p ,  o  ;  ]  }  , 

'  o  eo  i  =  l  °°  [kl'k2]  1 


k  1 '  2  ■ 


( V .  2 1 ) 


where  is  the  number  of  points  generated.  The  matrix  of  second 
partials  is  computed  using  standard  numerical  routines,  i.e., 
we  make  small  perturbations  in  dQ.  For  each  new  0  ,  we  recompute 
poles  G^,i=l,...,L  and  zeros  , i=l , . . . ,M,  which  occur  in  conjugate 
pairs.  Then,  we  obtain  the  corresponding  vectors  A  and  B  by  solving 

( 1-G^z  1)(1-G2z  1)...(1-Glz  =  1-a^z  1-...-aLz  L 


(V. 22] 


(1~P1Z  )(i-P2z  ^  *  *  *  ^ 1  ~ Pm2  ^  '  ' 


from  which  we  proceed  to  the  calculation  of  y(i)  and  eventually 


A[k  k  ]‘  Tke  suin  is  evaluated  for  the  several  perturba¬ 

tions  of  9q  needed  in  numerical  evaluation  of  the  matrix  of  second 

partials,  using  the  same  set  { C  C  r .  .  . ).}  each  time.  If  we  are 

1  9  ^  2  1 


interested  in  the  information  with  respect  only  to  poles,  we 

|G]_|  ,  .  .  .  , 

We  then  estimate  the  asymptotic  Cramer-Rao  bound  on  the  error 
covariance  matrix  for  0o  (note  that  we  consider  only  the  class  of 
asymptotically  unbiased  estimators)  as 


GL/2 '  '^1' ’  *  * ' ^L/2  ' 


use  the  reduced  parameter  vector  0 


a  a  m  /v  1 

E  <W  <VV  ;  I  <e0> 


(V. 23: 


where  the  matrix  ordering  is 

T  T 

R  >  Q  •*  y  Ry  >  y  Oy 


for  R  and  Q  square  of  dimension  (say)  K  and  v  any  vector  of 
dimension  K. 

Questions  regarding  the  details  of  the  Monte  Carlo  simulation 

are  under  investigation.  Next,  we  will  consider  the  special 

2 

case  of  AR  time  series  in  white  noise  of  variance  o  .  The 

n 

effect  of  this  noise  on  the  information  in  Cr,  .  is  evaluated 

l  K  ^  /  K  2  J 

n 

by  replacing  y(0)  obtained  from  (V.3a)  by  y’(0)  =  y(0)  +  . 

This  will  give  a  quantitative  measure  of  the  effect  of  the  model- 
change  phenomenon  in  the  presence  of  noise.  Studies  of  more  general 
time  series  and  noise  will  follow. 


58 


VI.  SPECTRAL  ESTIMATION  FOR  NOISY  COMPLEX  SINUSOIDS 

High  resolution  spectral  estimation  has  found  special 
popularity  in  applications  involving  real  or  complex  sinusoidal 
signals.  Examples  of  these  applications  include  radar  doppler 
processing,  and  radar  or  other  sensor  array  processing  for 
improved  angular  resolution  [27].  A  method  that  has  found 
special  appeal  in  these  applications  is  the  maximum  entropy 
method  of  spectral  analysis  (MEM)  introduced  by  Burg  [  8  ] . 

However,  there  have  been  several  disturbing  problems  with  this 
method  when  applied  to  sinusoids,  notably,  frequency  errors 
depending  on  the  sinusoidal  components'  phases  and  line  splitting 
under  high  order  estimates  and  large  signal- to-noise  ratios. 

We  have  considered  the  frequency  error  problem  and  have 
investigated  a  modification  to  Burg's  original  alaorithm, 
which  we  denote  the  tapered  Burg  algorithm.  The  tapered  Burq 
technique  is  a  direct  result  of  considering  a  weiqhted  least- 
squares  fit  to  the  parameters  of  the  all-pole  model,  subject 
to  Levinson's  recursion  constraint,  in  place  of  the  usual 
unweighted  least-squares  fit.  The  algorithmic  consequence 
of  this  approach  is  the  inclusion  of  an  appropriate  taper 
in  the  calculation  of  the  partial-correlation  coefficients  in 
the  usual  Burg  algorithm.  Based  on  the  expression  for  the 
frequency  error  in  the  spectral  estimate  of  a  sinusoid,  an 
optimum  taper  is  derived.  The  performance  of  this  optimum  taper 
is  then  compared  with  those  of  the  rectanqular  (untapered  Burg) 
and  Hamming  tapers.  It  appears  that  this  taper  makes  MEM  spectral 
estimates  of  sinusoids  using  Burg's  technique  more  robust,  without 
sacrificing  its  resolution. 
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In  the  popular  derivation  the  MEM  coefficients  define  a 
predication  error  filter  (PEF)  and  are  chosen  so  as  to  minimize 
the  average  of  the  forward  and  backward  prediction  residual 
energies  subject  to  Levinson's  recursion  constraint.  In  [8] 

Burg  generalized  this  derivation  by  minimizing  a  weighted 
average  of  the  residual  energies.  This  was  presumably  done  to 
allow  the  analyst  to  reflect  his  or  her  confidence  in  possibly 
disjoint  data  records.  For  contiguous  data,  it  did  not  appear 
that  any  weighting  of  the  average  residuals  was  needed. 

Therefore,  the  popular  MEM  spectral  estimate  has  been  that  of 
Burg's  original  suggestion  of  using  a  rectangular  (uniform) 
taper.  Recently,  Swingler  showed  [28],  through  numerical 
simulations,  that  a  reduction  in  the  error  in  the  estimated 
frequency  of  a  real  sinusoid  is  obtained  if  a  Hamming  taper  is 
employed  in  the  calculation  of  Burg’s  partial-correlation 
(PARCOR)  estimates.  This  method  of  tapering  is  exactly  that 
reported  in  [  8 ]  by  Burg. 

In  this  chapter,  we  first  derive  an  expression  for  the 
estimated  frequency  error  of  a  real  sinusoid  using  Burg's  tapered 
method.  This  is  a  generalization  of  the  error  expression  for  the 
un tapered  Burg's  technique  reported  by  Swingler  [29].  We 
subsequently  use  the  error  expression  to  derive  an  "optimum" 
taper.  Finally,  simulation  results  comparing  the  optimum, 

Hamming  and  rectangular  tapers  in  obtaining  Burg  spectral  estimates 
of  complex  sinusoids  in  noise  are  given. 

VI . 1  Generalized  Error  Expression 

In  this  section  we  derive  an  expression  for  frequency  error 
of  sinusoid  based  on  the  tapered  Burg  technique  for  a  aeneral 
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weighting  function.  Let  the  taper  wm<t)  be  defined  as  a 
function  of  the  continuous  variable  t  for  |t|  <  and  take 

on  the  value  zero  elsewhere.  Let  the  discrete  time  version  be 


,,  N-m. 

wmk  =  wm(k  "  "IT*  • 


(VI.l) 


As  in  (6)  we  assume  that  the  window  is  normalized  and  non-negative: 


N-m 

E  w  ,  =1  ,  all  m  , 

ink 


Wmk  ±  0 


,  all  m,k  . 


The  non-negativity  insures  that  the  magnitude  of  the  partial 
correlation  coefficients  (PARCORs)  do  not  exceed  unity. 

For  the  data  record  x.  ,k=0 , 1 , . . . , N ,  the  mth  PARCOR  is  given 
in  this  tapered  Burg  method  as  [  8  ] 


kt()  WmkDmkEmk 

a  =  -  2  rr - - - 

mm  N-m  .  ? 

k!0  ”mk(Dmk  *  Emk> 


(VI. 2) 


where  and  are  given  in  (VI.  6). 

For  x^  =  cos  (8k  +  <p)  ,  the  first  PARCOR  is 


a-^  =  - j cos (9)  + 


sin  (9)  E  w.,cos (20k+0+2$) 
k=0 
N-l 

l+cos(9)  E  w. ,  cos  ( 29k  +  6+2is ) 
k=0 


(VI. 3) 


We  choose  to  express  the  first  PARCOR  as 


=  -cos  (9  +  5) 

=  - [cos ( 0) +6sin  (0) ) 


for  5  <<  1  . 


(VI. 4) 
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From  (VI. 3)  we  identify  6  as 


N-l 

Z  w. ,  cos  ( 26k  +  0  +  2<}> ) 

k=0  iK _ 

N-l 

l+cos(0)  Z  w, ,  cos  (20k  +  0  +  2<j>) 
k=0 


(VI. 5) 


To  look  at  the  conditions  which  make  6  small  we  further  assume 
that  the  window  w^(t)  is  even  and  thus  its  Fourier  transform 
W^(oi)  is  real.  The  summation  in  (VI.  5)  then  becomes 

N-l 

Z  w, ,  cos  (2ek  +  9  +  24>)  = 
k=0 

(VI. 6) 

CO 

cos(N8  +  2i>)  Z  (-l)n(N_1)W1(2^n  +  2c!) 

n=-co 

Since  (o)  is  bandlimited  on  the  order  of  r/s  we  can  approximate 

the  summation  when  | -3  |  <  a  and  N  >>  1  as 


N-l 

Z  w. .  cos  ( 20k+9+24>)  -  cos  (N0  +  24>)  W.  (20) 
k=0  K  1 

-  0  f  ntx  <  ■ 9 1  <  11 


IT 


(VI. 7) 


From  (VI.  6)  and  the  constraints  on  w^  we  have  that  W^(0)  -  1. 
Thus  we  see  that  the  assumption  5  <<  1  is  valid  when  J  6  J  <  ~ 
and  N  >>  1. 

The  second  PARCOR  is  exactly  given  by 


a22 


A 

B 


1  + 


1  + 


B 

A 


N-2 

Z 

k=0 


A 

B 


N-2 

Z 

k=0 


W2j{cos  (2  0k  +  29  +  2i}>) 


W2^cos(20k  +  20  +  2.t>) 

V 


(VI. 8) 


i 

i 


where 
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2 

A  =  a^  +  2a^^cos(6)  +  cos(29) 

2 

B  =  +  2a^cos{0)  +  1 

Using  the  approximation  for  results  in 

A  =  -  ( 1- 6 2 )  sin  (6)  ,  B  =  ( 1- 5 2 )  sin  ( 0 )  , 

and  thus 

a 22  =  1  »  5  <<  1  . 

The  Levinson  recursion  gives  a2^  =  a^-,  +  a22all  ~  '2all  an<^  t*le 
resulting  PEF  is  then  (1,  -2  cos  (0-6)  ,  1).  Comparing  this  to 
the  ideal  PEF  (1,  -2  cos  (0)  ,  1) ,  we  identify  the  frequency 
error  as 


N-i 

E  w,,  cos  (29k  +  e+24>) 

1  k-ft 

=  sin  (  0 )  — - j -  (VI.  9) 

l+cos(9)  E  w. ,  cos  ( 2  6k  +  9  +  2<J> ) 
k=0  iK 


which  may  be  approximated  by 


if 


For 


N-l 

E 

k  =  0 


1 

-  —  sxn(9)  E  w,  ,cos  ( 20k  +  0  +  2.f ) 

k  =  0 

FFT  <  1 6  |  <  tt  - 

-  sin (0) cos (NO+2^)  (29)  . 

uniform  weighting,  i  ,  we 

wlkc°s  (20k  +  0+2<{>)  =  i  cos(N0  +  2if>) 


have 

sin (N9 ) 
sin  (0) 


(VI.10 (a) ) 

(VI. 10 (b) ) 


and  (10a)  reduces  to  Swingler ' s  expression  [29]. 
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VI. 2  An  Optimum  Taper 

If  we  consider  the  phase  ^  as  a  random  variable  uniformly 
distributed  on  [  —  rr ,  tt  )  then  the  mean  value  of  the  frequency  error 
is  zero,  using  (10a) .  The  variance  of  this  frequency  error  is  then 


,  ?  N-l  N-l 

var(Af)  =  — =■  sin  (0)  Z  Z  w, .  w.  .cos  (29  (k-  i)  ) 
8:r2  k  =  0  £=0  ik  11 


(VI. 11) 


As  a  criterion  for  selecting  a  taper  we  use  the  average  frequency 


error  variance: 


<var(Af)>  =  —  !  var(Af)d6 


N-l  N-l 


L  W,,W 


5,.  .  .  ~  r  5 


8  tt  k  =  0  1=0 


lk"ll[2  °k-£  4  ul-k+l  4  ul+k-!j  ' 


(VI. 12) 


where  the  subscripted  delta  is  the  digital  (Kronecker)  impulse. 

With  the  normalization  constraint  introduced  using  the  Laarange 
multiplier  X  the  resulting  optimum  taper  is  given  by  the  solution  of 


-  -1  -  '  Cii  2 ,Ci j 


1  ,  j  =  iil 
0  ,  otherwise 


(vi.  i3: 


The  general  system  of  equations  with  a  tri-diagonal  coefficient 
matrix  has  a  known  recursive  solution  related  to  the  LU  decomposition 
of  the  coefficient  matrix  (see  for  example  [30]).  For  the  system 
of  equations  with  the  special  coefficient  matrix  in  (VI. 13) ,  we 
have  derived  closed  form  expressions  for  the  taper  coefficients 
(Appendix  B) .  These  are  given  for  the  m-th  order  PARCOR  by: 


6  (k-t-l)  (N-m-k  +  1 ) _ 

;mk  (N-m+1) (N-m+2) (N-m+3) 


,  k=0 , . . . , N-m 


(VI .  14) 
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Thus,  this  optimum  taper  is  parabolic  in  form,  it  is  even, 
positive  and  has  a  maximum  at  k  =  N  ~  •  Furthermore,  knowinq 
wmQ  and  w  ^  from  (VI.  14),  one  can  generate  the  remaining 
coefficients  recursively  from: 


w  .  —  2w  ,,  —  w  , ,  ..  —  \ 

mk  m(k-l)  m(k-2) 


(VI. 15) 


VI . 3  Simulation  Results 

The  effect  of  tapering  was  numerically  investigated  by 
generating  spatial  samples  of  complex  sinusoids  in  noise  for 
various,  phase  and  angle  (frequency)  combinations.  For  brevity 
of  space,  however,  only  the  results  of  a  very  few  representative 
simulations  are  included  in  this  report.  The  method  designations 
on  the  plots  are:  Burg  for  untapered  Burg  method,  WBurgH  and 
WBurgO  for  tapered  Burg  technique  using  the  Hamming  taper  and 
the  optimum  taper  respectively.  The  Hamminq  taper  is  used  for 
comparison  with  the  simulations  in  [28] . 

Figure  (VI. 1)  shows  three  tapered  Burg  spectral  estimates 

of  two  complex  sinusoids  at  -30°  and  -45°  off  broadside  with  zero 

phases  each,  in  complex  white  noise.  The  signal- to-noise  ratio 

is  15  dB  where  SNR(dB)  =  10  log  )  .  it  is  obvious  in  this 

on 

example  that  WBurgO  has  the  smallest  frequency  errors  as  well  as 
the  highest  resolution.  This  has  been  the  case  in  the  great 
majority  of  the  examples  that  we  have  run  thus  far.  Except 
for  a  few  cases  the  bias  of  WBurg  optimum  has  been  lower  than  Burq. 

The  resolution  of  WBurg  Hamming,  however,  was  found  to  be  consistently 
poorer  than  Burg  and  WBurgO.  Figure  (VI. 2)  shows  the  spectral 
estimate  of  a  real  sinusoid  with  a  large  SNR.  It  can  be  seen  that 
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WBurg  optimum  taper  has  reduced  the  frequency  error  and  prevented 
peak  splitting  which  is  shown  by  the  untapered  Burg  method  of 
order  4.  This  again  was  accomplished  without  sacrificing 
resolution . 

Figure  (VI. 3)  is  a  plot  of  the  maximum  frequency  error  for 

each  method  for  a  unit  amplitude  real  sinusoid  at  angle  (frequency) 

of  0  =  t/5  as  a  function  of  the  number  of  data  samples.  The 

maximum  frequency  errors  were  found  by  solvinq  for  the  pole 

locations  over  the  range  of  values  for  the  sinusoid  phase  <J> . 

2 

The  noise  variance  is  a  =  0.025.  It  should  be  noted  that  the 

n 

optimum  taper  was  derived  on  the  basis  of  minimum  average  frequency 
variance .  This  simulation  shows  comparative  error  reduction  for 
the  maximum  frequency  error  only  for  0  =  rr/ 5  between  the  optimum 
and  Hamming  tapers.  As  was  pointed  out  earlier,  however,  all 
of  our  simulations  have  shown  better  resolution  for  the  optimum 
taper  compared  to  the  Hamming  one.  It  is  also  interesting  to  note 
that  the  main  reduction  in  frequency  error  is  apparently  obtained 
by  tapering  the  first  order  residual  enerqy  only.  This  has  also 
been  borne  out  by  the  majority  of  our  simulations. 


xponentials 

inusoid)  in 
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VII.  SUMMARY  AND  CONCLUSIONS 

Data-adaptive  spectral  estimation  methods  have  found 
widespread  application  in  recent  years,  notably  in  radar  signal 
processing.  The  work  summarized  in  this  report  was  concerned 
with  two  such  methods:  the  autoregressive  (AR)  and  the  auto¬ 
regressive-moving  average  (ARMA)  spectral  estimators.  Our 
main  concern  in  this  investigation  was  the  resolving  capabilities 
of  the  different  models  in  the  presence  of  noise  and  the 
statistical  properties  of  some  new  spectral  estimators. 

Since  the  AR  methods  are  by  far  the  most  popular,  we 
postulated  an  AR  noiseless  signal  model.  The  resulting  appropriate 
model  in  the  presence  of  noise  was  then  shown  to  be  ARMA.  To 
investigate  the  degradation  in  resolution  in  the  spectrum  of 
noisy  signals,  we  also  used  a  parallel  resonator  model  for  the 
signal.  This  model  was  also  shown  to  be  equivalent  to  an  ARMA 
one.  The  resonator  model,  however,  made  it  possible  to  evaluate 
the  loss  in  resolution  as  a  function  of  signal- to-noise  ratio  and 
signal  parameters. 

Since  an  ARMA  model  may  be  approximated  by  a  high-order 
AR  one,  attention  was  next  focussed  on  the  development  of  a 
robust  method  for  identifying  the  order  of  an  AR  model.  A  method 
closely  related  to  that  of  Akaike ' s  information  criterion  was 
developed.  This  method  proved  to  be  more  stable  than  the 
minimum  AIC. 

ARM^i  spectral  estimation  was  treated  extensively  in  this  work. 
The  generality  of  the  model  was  important  from  two  distinct  points 
of  view:  i)  as  the  appropriate  model  for  noisy  AR  signals,  thus 
to  improve  resolution;  ii)  as  a  resonable  model  for  signals  with 
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spectra,  narrow  or  wide,  possessing  deep  nulls  and/or  sharp 
roll-offs.  Because  of  the  computational  complexities  of  the 
optimum  (maximum  likelihood)  method  of  estimating  the  ARMA 
parameters,  several  sub-optimum,  computationally  efficient 
techniques  were  devised.  These  methods  behaved  reasonably  well 
for  moderate  to  large  data  samples,  but  were  inferior  to  the 
optimum  one,  as  expected,  in  terms  of  statistical  efficiency. 
The  efficiency  of  various  classes  of  ARMA  spectral  estimators 
are  currently  under  further  investigation  by  Monte  Carlo 
simulations  and  will  be  reported  on  in  the  future. 

The  maximum  entropy  method  (MEM)  of  spectral  estimation, 
using  Burg's  technique,  has  been  very  popular  in  estimating 
the  spectra  of  sinusoidal  type  signals.  It  has  been  known, 
however,  that  the  accuracy  of  such  spectral  estimates  is 
significantly  influenced  by  the  phase  of  the  sinusoids.  Based 
on  some  recent  results,  we  developed  an  optimum  taper  for  the 
residual  enerciies  in  the  Bure  MEM  technique.  It  was  shown  that 
the  use  of  such  a  taper  substantially  reduced  the  sensitivity 
of  the  spectral  peaks  to  the  sinusoids'  phases  and  markedly 
reduced  the  occurrence  of  line-splitting  in  the  usual  MEM 
estimates.  It  is  remarkable  that  these  improvements  were  made 
without  sacrifice  in  the  resolution  of  the  spectral  estimates. 
Work  is  continuing  to  more  fully  explain  the  effect  of  such  a 
taper  and  its  ramifications  in  Burg-type  parameter  estimation 
techniques . 


ni:  i 
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APPENDIX  A 


Equation  (11.25)  is 


r  (0)r  (2)  -  r;  (1) 

1  m _ Mn _ 1  m 

r2  (0)  -  r2  (1) 

V  V 

J  m  2  m 


(A .  1 ) 


Now  substitute  for  r  (0),  r  (2)  and  r  (1)  in  A.l, 


(aL°uacosV°w'  }  (amo°ua4cos  (2wnf  V  >  "  (anoa-aaDmcos  (wnf  V  5  2 


2  2  2  2  2  2  2 
(a  a  acosii  +o  ,)  -(a  o  a,o  cos  (w  -4  )) 
mo  u  Tm  w  mo  u  m  m  vm  ' 


(A. 2) 


Expanding  the  numerator  and  denominator  of  equation  A. 2,  we  aet, 


22  22  2222 
(a  c  )  a  P  cospcos  (2w  -4  )  +a  o  a  ,  a r,  cos(2w  -4  ) 
mo  u  m  m  rm'  mo  u  w  m  m  m 


2  2  2  2  2 
-  (a  o  )  a  o  cos  (w  -4  ) 
mo  u  m  m  m 

2 

22  2  2  4  222 

(a  a)  a  cos“4  +cr  ,+2a  a  a  ,acos4 
mo  u  m  w’  mo  u  w’  vm 


2  2  2  2  2 
-  (a  o':)  ct  o  cos  tw  -4  ) 
mo  n  m  m  m 


(A. 3) 


Define  the  signal- to-noise  ratio  ?  as  , 


ry  (0>  "  °w- 
1  m 


r  (0) 
m 


where , 


=  Y  -  1 


2  2 

r  (0)  =  a  a  ctcos4> 
y  mo  u  m 

m 


Dividing  the  numerator  and  denominator  of  A. 3  by  o^,  and  using 


(A. 4) 


A.  4,  we  get, 


'  2  cos(2wm-<,[n)  cos(2„m-»m)  2  2 

Y  COS«n>  '  COs4V  ‘  1  cos2f 

^  m 


2  2, 

-  PCOS  (W  -i  ) 

Y2  +  1  +  2y  -  y2  - -H-HL. 

cos^t 

m 


Rearranging  terms  in  A. 5,  we  have 


Y2  (cos$incos  (2wra-$m)-cos2  (wm-0m)  ]  +Ycos4>mcos  (2wm-0m) 

2  2  2  2  2 
Y  cos  <p  -p  cos  (w  -ip  )  +2ycos  <p  +c os  <t> 
m  ^m  m  m  1  m  m 


(A. 5) 


(A. 6) 


As  y 


->  X) 


1  and  the  left-hand  side  of  equation  A. 6  is, 


(14  ] 


cos 4>cos  ( 2w- 4j  ) -cos  (w  -<b  ) 
m  m  m  m  m 

2  2  2 

cos  4>  -o  cos  (w  —  <fe _ ) 

m  Mm  m  m' 


=  -1 


2  2 

Solving  for  p  cos  (w  - <b  )  in  A. 7  we  get, 
mm  m  ^ 


Substitute  A. 8  in  A. 6  to  get, 


Y2  cos^cos  (2wm-$m)  -cos2  (w^J  ]  +Ycos4»mcos  (2^-^) 

2  2  I  2  2 

Y  -cos^cos  ( 2wm—  <J>m)  +cos  lw  -4>m)  +2ycos  4>m+cos  ij> 


But  the  signal- to-noise  ratio  is  (equation  A. 4) 


5  =  Y-l 


A.  7  then  becomes. 


(£+1)  cos^cos  (2wm-4>m) -cos2  +(C+l)cos4>mcos  (2wm-4)m) 


(5  +  D  2  |  — cosmos  (2wm-.j >m)  +cos 2  |  +2  ( 5  +  1)  cos2<(>m+cos2<t>m 


(A. 7) 


(A. 8) 


(A. 9) 


(A. 10) 


Equation  A. 10  simplifies  to. 


74 


APPENDIX  B 

In  this  appendix  we  derive  a  closed  form  expression  for 
the  solution  of  the  following  equations: 


2  -1  0  ...  0 

W1 

'a' 

-1  2-1  0  ...  0 

• 

0  -1  2-1  ...  0 

= 

0  ...  0  -1  2 

A 

l  ) 

NJ 

(B.l) 


where  A  is  the  normalization  factor  for  W. 

A  recursive  algorithm  for  a  set  of  linear  equations  with 
a  general  tridiagonal  coefficient  matrix  is  given  in  [30].  We 
use  the  notation  used  in  this  reference  to  derive  our  closed 
form  solution. 

The  solution  for  w^  =  wN  is  given  by 
^N 


wM  =  w.  = 

N  la 


(B.2) 


N 


where  according  to  [30 


,  (k-1)  .  , 

gk  —  k  gk_ ^  ^  1/  ...  , N 


ak  =  2+ek  '  k=2. 


=  2 


•  ,N 


( B  .  3 ) 


6k  a 


k-1 


We  first  derive  an  expression  for  6^.  Let  8^ 
and  6^  are  integers. 

From  (B. 3) 

1 


^k  = 


2  +  8 


k-1 


=  -vk/6k/  where 


,  k  =  2 ,  ...  , N 


(b.  4: 
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Then 


\-l  ~  Vk- 1 


( B  .  5 ) 


Now  let  then 


fik-l  +  mk-l 


(B.  6 ) 


But  we  have  v2  =  1,  S-,  =  2  and  m2  =  ^  anc^  from  (B.6)  it  is 
obvious  that  =  constant.  Therefore, 


k- 1+1 


or  with  v2  =  1  ,  52  =  2 


(k—  1 ) 


(B.7) 


now  simply  follows  as 


“k  *  2tsk  -  ^  ■  k-1- 


(B .  8 ) 


Let  =  kg^  and  substitute  for  in  (B.3).  h^  then  satisfies 


the  difference  equation 


hk  "  hk-l  =  kX  '  ho  =  0 


(B .  9 ) 


Z-transform  (B.9)  to  get 


H  ( z )  = 


Z_1  (z-1)2 


(B. 10) 


resulting  in 
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h,  =  (step)  *  (ramp)  A  =  A  T.  j 

R  j=l 

_  k  (k+1)  , 


(B. 11) 


Therefore 


gk  2 


( B  .  1 2  ) 


Substitute  in  (B.2)  to  get 


NA 

WN  "  W1  ~  ~T 


(B  .  1  3 ) 


Using  (B.13)  as  the  initial  condition  we  can  now  solve  for 
other  w. 's,  noting  that 


Wk  “  2  wk-l  +  wk- 2 


=  -A  ,  k=2,  ...  , N 


(B . 14 ) 


wo  =  0  '  W1 


Z-transform  of  (B.14)  results  in 


w(z)  =  +  . 

(z-l)J  2(Z-1)Z  (2-1)  ^ 


(B . 15 ) 


Taking  the  inverse  transform  of  (B.15)  gives 


Lih  +  Ugm.  +  +  (k+1)  A 


..  _  k  (N+l-k )  , 

wk - 2 -  X 


(b. i6: 


The  normalization  factor  A  can  now  be  found  as 


A  =  (  £  Wt  ^ 
k=l  * 


(B . 17 ) 


Alternatively, 


A  =  \  . 

N+l  N  N+3 


(B. 18) 
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